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Summary. We study the formation and the dynamics of correlations in the ve- 
locity field for ID and 2D cooling granular gases with the assumption of negligible 
density fluctuations ("Homogeneous Velocity-correlated Cooling State", HVCS). It 
is shown that the predictions of mean field models fail when velocity fluctuations be- 
come important. The study of correlations is done by means of molecular dynamics 
and introducing an Inelastic Lattice Maxwell Models. This lattice model is able to 
reproduce all the properties of the Homogeneous Cooling State and several features 
of the HVCS. Moreover it allows very precise measurements of structure functions 
and other crucial statistical indicators. The study suggests that both the ID and 
the 2D dynamics of the velocity field are compatible with a diffusive dynamics at 
large scale with a more complex behavior at small scale. In 2D the issue of scale 
separation, which is of interest in the context of kinetic theories, is addressed. 

1 Introduction 

The inelastic hard spheres model [1] without energy input, initially prepared 
in a homogeneous state, exhibits, after a rapid transient, a regime charac- 
terized by homogeneous density and a probability distribution of velocities 
that depends on time only through the total kinetic energy (global granular 
temperature T g (t)), i.e. a scaling velocity distribution. This is the so called 
Homogeneous Cooling State (HCS) or Haff regime [2]. It has been shown 
by several authors [3-6] that this state is unstable with respect to shear and 
clustering instabilities: structures can form that seem to minimize dissipation, 
mainly in the form of velocity vortices and high density clusters. These insta- 
bilities grow on different space and time scales, so that one can investigate 
them separately. Several theories have been proposed to take into account 
the emergence of structures in granular gases. Some of these are more fun- 
damental, because derive the correlation functions directly from the kinetic 
equations [7]; others, that assume the validity of the hydrodynamic descrip- 
tion, deserve the name of mesoscopic theories [8] ; others are phenomenolog- 
ical theories that suggest analogies with Burgers equation [9,10] or spinodal 
decomposition models [11] or mode coupling theories [12]. Some of these the- 
ories can describe the behavior of the cooling granular gas far deeply into the 
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correlated regime, giving predictions for the asymptotic decay of energy. In 
this paper, after a brief summary of the results of these theories, we show 
how the departure from the HCS can be well modeled by a class of mod- 
els obtained placing on a lattice the original homogeneous Inelastic Maxwell 
Model. These models have the disadvantages of being conceived under the as- 
sumption of negligible density fluctuations and therefore can be useful only in 
the description of the first instability i.e. the growth of velocity fluctuations. 
In section 2 we review the HCS, its instabilities and the existing theories. In 
section 3 we briefly discuss the Inelastic Maxwell Model [13-15], which is a 
starting point for the introduction of the Inelastic Lattice Maxwell Models 
(ILMM). In section 4 and in section 5 the analysis and results of the ILMM in 
one and two dimensions are reviewed [14,16]. Finally in section 6 conclusions 
are drawn. 

2 Instabilities of the Homogeneous Cooling State 

A cooling granular gas in d dimensions is defined as an ensemble of N grains, 
i.e. hard objects (rods if d = 1, disks if d = 2, spheres if d = 3) of linear size 
(diameter) a, placed in a volume V with periodic boundary conditions. The 
grains evolve freely and interact with each other through instantaneous binary 
inelastic collisions. The rule that gives the velocities after the collision as 
functions of the velocities before the collision is the definition of the particular 
granular gas model. In this case we use the model with constant restitution 
coefficient without rotational degrees of freedom and set the mass of the 
grains m = 1. Other models have been considered in the literature [1]. The 
collision rule between a particle with velocity v and one with velocity v* for 
this model is: 



where the primed velocity is the post-collisional one, n = (r — r*)/|r — r*| is 
the unit vector in the direction joining the centers of colliding particles, and 
r is the restitution coefficient and takes values between and 1. When r = 1 
the collision is elastic. 

Usually (in numerical or real experiments) granular gases are prepared in 
a homogeneous situation: uniform-random positions of grains, Gaussian or 
uniform-random initial velocity with no preferred direction. It happens that, 
if the system is large enough or the inelasticity large enough, the imposed 
homogeneity is broken after a certain time. The more accepted scenario is 
a two time homogeneity breaking: at a time t s the velocity field becomes 
unstable to the formation of shear bands, then at a time t c > t s the density 
field becomes unstable toward the formation of high density clusters. 




(1) 
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2.1 The homogeneous cooling state 

A granular gas evolving from a homogeneously random state loses memory of 
its initial condition after a time of the order of one collision per particle and 
rapidly enters the Homogeneous Cooling Regime. This regime is expected to 
be well described by the granular Boltzmann Equation (see for example [7]). 
This means that in this regime the total probability density function (p.d.f.) 
at collision is factorisable in one-particle p.d.f.'s: pjv( r i-- r iv, Vi...vjv, t) = 

The kinetic definition of HCS is given by the homogeneity ansatz plus the 
scaling ansatz for the one particle distribution function: 



p l(r ,v,t) = -P(v,t) = — — P(c) (2) 



where c = v/vo(t) and vo(i) is the thermal velocity defined by T(t) = Vo(t)/d 
with T(t) the granular temperature; here we have assumed that J dv J drpi = 
1 and also that J dvP = 1. If the Eq. (2) is inserted in the Boltzmann 
Equation an equation for the temperature is obtained: 

f - " 2 ^ T < 3 > 
where u> is the time dependent collision frequency, while 7 is the time 
independent cooling rate. These two functions can be approximated, using 
the Maxwellian approximation P w (^)~ d / 2 exp(— ^f-), by loq and 70: 

loo oc VT (4a) 
1 -r 2 

70 = — • (4b) 
In this case the solution of the temperature equation (3) reads: 



nt) = r, , :!u v - cxp(-2 7o t) (5) 



T(0) 
(1 + 7 ot/to) 2 

where to = l/w(T(0)) is the mean free time at the initial temperature 
T(0) and 

r = — ln(l + 70 + t/t ) (6) 

7o 

is the cumulated collision number obtained from the definition dr = u)(T(t))dt. 
Eq. (5) is known as Haff 's law [2] . 

Corrections to the constants appearing in eq. (5) stem from a more care- 
ful consideration of the HCS: when the volume fraction is non negligible the 
Enskog-Boltzmann equation should be employed instead of Boltzmann equa- 
tion. This is identical to the Boltzmann equation but for a multiplicative 
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Fig. 1. Growth rates Ca/7o for shear (A =_L), heat (A = H) and sound (A = ±) 
modes versus ka for inelastic hard disks with r = 0.9 at a packing fraction cf> = 0.4. 
The dashed line indicate the imaginary parts of the sound modes. Here a is the 
diameter of a particle. (From Orza et al. [19]) 



constant in the collision integral that takes into account static density cor- 
relations due to the fact that the gas is not perfectly dilute. Very recently 
it has been shown [17] that there are also small velocity correlations which 
must be considered and that modify the molecular chaos hypothesis: with 
this correlations the law T ~ t~ 2 is still valid but more precise corrections to 
the constant 7o/io appear. 

The Haff law can be also derived in the framework of granular hydro- 
dynamics [18,6]: in the HCS hydrodynamics can be considered valid as a 
consequence of homogeneity of the density and velocity field and of the slow 
temperature decay. The evolution of the (space uniform) temperature field is 
exactly the same as in eq. (5). 

2.2 Instabilities of the homogeneous cooling state 

Kinetic arguments, linear stability analysis of hydrodynamic equations and 
numerical simulations lead to the evidence that the HCS is unstable. In par- 
ticular, it has been shown that large scale fluctuations of velocity and density 
can grow exponentially breaking the homogeneity of the system. The fact that 
the fluctuations arise at large scale means that, if the system is small enough, 
the HCS lasts forever. Moreover the minimum scale of formation of fluctu- 
ations depends upon inelasticities and is smaller at larger inelasticities: this 
means that for a quasi-elastic gas it is necessary to take a very large volume 
to observe the breaking of homogeneity. 
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The linear stability analysis of hydrodynamics gives a good description of 
the instability of fluctuations. The procedure is applied to the fields rescaled 
with respect to the Homogeneous Cooling State, i.e. n = n(r,t)/n (being n 
the average density), u(r,i) = u(r, t)/y/T{t) and f(r,t) = T(r,t)/T(t), and 
therefore any instability or emerging structure is not absolute but relative to 
this decaying reference state. Moreover the instabilities are studied in k space 
(the space of Fourier modes) and, when d > 1 , the velocity field is decom- 
posed in a parallel component uy (k, t) and in d — 1 orthogonal components 
Uj_(k,t) with respect to the vector k. In figure 1 we report the results of 
such an analysis, taken from the literature [19]. This analysis shows that the 
evolution of fluctuations of normal velocity components (shear modes, uj_) 
are not coupled with any other fluctuating component. At the same time, 
the remaining components are coupled together. The rate of decay/growth 
for the shear mode reads C-l(^) = 7o(l — ^ 2 £i_) where £j_ depends on the 
transport coefficients appearing in the hydrodynamics. We refer to [8] for 
detailed calculations of this correlation length. At low values of k (in the so- 
called "dissipative range" ) also the heat mode is "pure" , as it is given by the 
longitudinal velocity mode uy only, with eigenvalue Cn(k) ~ 7o(l — £ 2 fc 2 ); in 
this range the sound modes are combination of density and temperature fluc- 
tuations. The most important result of this analysis is that Ci(&) and Cff(&) 
are positive below the threshold values kj_ — l/£j_ ~ \ft and k* H ~ 1/£|| <~ e 
respectively, indicating two linearly unstable modes with exponential (in t, 
i.e. the time measured by the cumulated number of collisions per particle) 
growth rates. Here we have used the notation e = 1 — r 2 = 2djo. 

The shear and heat instabilities are well separated at low inelasticity (i.e. 
low e), as k*^ ~ sfi while k* H ~ e, so that kj_ 3> k* H . It is also important 
to note that the linear total size L of the system can suppress the various 
instability, because the minimum wave number k m i n — 2tt/L can be larger 
than k* H or even than k*\_. 

When fluctuations grow, structures emerge. In molecular dynamics sim- 
ulations these structure appear as shocks in d = 1, or vortices and clusters 
in d = 2. In > 2, a detailed analysis of structure factors (from fluctuat- 
ing hydrodynamics [8] or ring kinetic theory [7]) has established that all the 
main correlation lengths (which indicate the typical size of structures in the 
rescaled fields) should grow as ~ t 1 / 2 . In particular this is expected for the 
size of vortices L v <~ t 1 / 2 and the size of clusters L c i <~ (r/e) 1 / 2 : this again 
shows that, at small inelasticities e, the size of clusters grows more slowly 
than that of vortices. 

2.3 Conjectures 

After a brief transient during which the instabilities grow exponentially in 
r with the rates given above, the description suggested by linear stability 
analysis ceases to be meaningful. Molecular dynamics (MD), in the absence of 
experiments (the cooling granular gas being a reference model and not a real 
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system), are the more reliable tools, but they are very cpu-time consuming: 
to obtain a decent statistics of the asymptotic dynamics of a granular gas, 
large enough to observe instabilities, MD simulations must be carried on for 
days. To our knowledge a definitive clear picture of the asymptotics of this 
model is still lacking (with month-lasting simulations Isobe [20] has obtained 
what he called the "final state" , which resembles 2D turbulence, but this was 
done in a quasi-elastic limit). 

On the other side there are effective models that can grasp the essence of 
the phenomena with particular assumptions. Here we rapidly review some of 
these: 

• Brito and Ernst have shown by means of a mode-coupling [12] theory 
that the asymptotic energy decay should follow a diffusive form ~ T~ d l 2 . 

• Wakou et al. [11] have demonstrated that the evolution of the flow field 
of a granular fluid, neglecting the convective term u • Vu, can be cast 
in the form of a Time Dependent Ginzburg-Landau equation for a non- 
conserved order parameter. The energy functional has a continuous set 
of degenerate minima, having the shape of a Mexican hat. When the 
convective term is added, only subsets of admissible solutions are selected 
out of this infinite set of minima. In two dimensions only two distinct 
minima survive: therefore the d — 2 cooling granular gas, during the 
formation of instabilities, greatly resembles the spinodal decomposition 
for a non-conserved order parameter (Model A universality class, in the 
Hohcnbcrg-Halperin classification [21]). 

• Ben-Naim et al. [9] have studied the cooling granular gas in d = 1, dis- 
covering that, after the homogeneous (Haff) phase, it asymptotically be- 
comes independent of the value of inelasticity e and maps onto the sticky 
gas model. For such a sticky gas the asymptotic temperature decay is pro- 
portional to i~ 2 / 3 , when density clusters and velocity shocks form. The 
velocity field of the sticky gas [22] is described by the inviscid limit of the 
Burgers equation [23]. The relation to the Burgers equation is useful to 
establish an estimate of the tails of the asymptotic velocity distribution 
P(v, t) ~ ^/^(ut 1 / 3 ) (independent of r), which reads $(z) ~ exp(-|z| 3 ). 
The MD simulations have revealed only slight deviations from the Gaus- 
sian, and the authors have imputed the discrepancy from the expected 
behavior to the smallness of the constant in front of |z| 3 . More recently 
it has been conjectured [10] that the Burgers equation describes the flow 
velocity u(r, t) of a cooling granular gas for arbitrary values of the inelas- 
ticity e in generic d dimensions. This conjecture implies that the asymp- 
totic behavior of the cooling inelastic gas is independent of e (as it always 
falls in the universality class of the sticky gas) and that the upper crit- 
ical dimension for the disappearance of the inelastic collapse is d c = 4. 
However MD simulations in d > 4 seem to show that the shear instability 
and the inelastic collapse do not disappear [24]. 



Cooling granular gases 7 



3 A starting point: the homogeneous Inelastic Maxwell 
Model 

To better illustrate the analysis of the Inelastic Lattice Maxwell Models 
(ILMM) we start with a brief description of the homogeneous Inelastic Maxwell 
Model (IMM). The definition of the model is given in terms of its Boltzmann 
Equation, which, for scalar velocities (d = 1), reads: 



<9t 



(7) 



with 7* = (I — r)/2. The IMM can be introduced in different ways: (a) 
its Boltzmann equation can be regarded as the master equation of a simple 
stochastic process inspired to a model put forward by S.Ulam [25]: the pro- 
cess consists in the evolution of an ensemble of N velocities, at each time 
step (corresponding to At = 2/N, so that r increases of 1 after averagely 
one collision per particle) two velocities are chosen at random with uniform 
distribution and are changed according to the inelastic collision rule; (b) the 
associated Boltzmann equation can be considered the inelastic generaliza- 
tion of the equation for elastic Maxwell Molecules [26,27]; (c) its Boltzmann 
equation can be derived from the ordinary Boltzmann equation for inelastic 
hard spheres assuming that the term \v — v'\ in the collision integral can be 
approximated by VT, as suggested by Bobylev et al. [28]. 

The recent surge of interest for this model has been triggered by the 
discovery of an exact solution [13-15] for the d = 1 case (eq. (7)) and the ap- 
pearance of computable power law tails in the solutions of d > 1 cases. [29,30]. 
In a previous paper Ben-Nairn and Krapivsky [31] showed that eq. (7) could 
not have a scaling solution with finite moments: in fact, a scaling form 

PKr)--|r/(-f-) (8) 
vo{t) Vo(T) 

with Vq(t) = J v 2 P(v,r)dv — E{t) ~ cxp(— ra 2 ), imposes a temporal de- 
pendence of the higher moments (if finite): 



(v 2m )=v ( T ) 2m ^ m (9a) 

where M2m = J dyf(y)y 2m (9b) 

does not depend on time. But they found that 

(v 2m (T)) 

lim n 2, \ = °°- ( 10 ) 

This result does not exclude the existence of a scaling solution (8) but 

requires a solution with diverging moments of order 2m > 4, i.e. very large 
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Fig. 2. Asymptotic velocity distributions P(v,t) versus v/vo(t) for different values 
of r from the simulation of the Inelastic Maxwell Model in ID (left) and 2D (right). 
In ID the asymptotic distribution is independent of r and collapse to the Eq. (11). 
In ID, the chosen initial distribution (exponential) is drawn (same result with 
uniform and Gaussian initial distribution). In 2D the distributions still present 
power-law tails, but the power depends upon r: for r — > 1 the pdf tends to a 
Maxwell distribution. Data refers to more than N = 10 6 particles. 



algebraic tails. In this case standard methods [27] based on Fourier transform 
fails due to the singularity in k = corresponding to such algebraic tails. 

Remarkably, a direct inspection shows that the following velocity distri- 
bution 



P(v,t) = 



ttv {t) 



(11) 



is indeed a solution of the non-linear Boltzmann equation (7), for every value 
of r [13-15]. The solution (11) has the property that the moments of order 
2m > 4 diverge, and its Fourier transform has a singularity of the type \k\k 2 
ink = 0. 

We believe that cq. (11) represents the asymptotic solution for a large 
class of initial distributions, for the following arguments: 

• as shown by Ben Nairn and Krapivsky [31], the dynamics of the moments 
for a generic starting distribution can be computed, giving the following 
limit lim t ^ 0O (w 2m (t))/(w 2 (t)) m = oo for m > 1; 

• secondly we performed numerical simulations of the BK model, collecting 
evidence of the convergence to the solution (11) for several starting veloc- 
ity distributions, namely uniform, exponential (see figure 2, left frame) 
or Gaussian. 
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For higher dimensions, the problem of a scaling solution has been recently 
addressed in [10,29]. Assuming a solution with large algebraic tails, it is pos- 
sible to compute the asymptotic exponent as a solution of a transcendental 
equation. Such equation depends explicitly on the restitution coefficient, indi- 
cating that the tail exponent varies with r (and correctly diverges for r — > 1 , 
corresponding to the recovering of the Maxwell distribution for the elastic 
case). 

A numerical observation of such algebraic tails in 2 dimensions are shown 
in right frame of figure 2, where a dependence on r is put in evidence. A more 
precise comparison of the measured exponent and the theoretical predictions 
are shown in figure 3. In this figure we show the results of extensive numerical 
simulations for the completely inelastic case (r = 0) of several models, in two 
and three dimensions. 

In fact, as noted in [29], the tail exponent is in general dependent on 
the details of the collisions, as for instance on the distribution of the impact 
parameter. Two different choices have been studied in 3 dimensions, corre- 
sponding to the models denoted IMM- A and IMM-B in reference [29] : in the 
first model grazing collisions are favored with respect to the second one. It is 
interesting to notice that the analytical and numerical study confirm the ex- 
istence of a scaling solution with algebraic tails, irrespectively of the collision 
details, but that the precise value of the exponent is not universal. 

The relevance of such a scaling solution for d > 1 has been addressed by 
Bobylcv et al. [32], who have proved the existence and the uniqueness of such 
an asymptotic solution (giving details on the basin of attraction). 

Interesting features emerge when the inelastic Maxwell model is extended 
to treat grains having different physical properties, such as unequal masses, 
different restitution coefficients, radii, etc. The study of a binary Maxwell 
mixture with scalar velocities showed few unexpected properties [33,34]: 

1. the kinetic temperatures of the two species are different, i.e. there is no 
energy equipartition, but the the ratio Ti/T 2 reaches asymptotically a 
constant value. Such a feature agrees with the results obtained in the 
framework of the Boltzmann-Enskog transport equation by Garzo and 
Dufty [35]. 

2. The velocity distributions in general differ in shape. 

3. The velocity distributions have power law tails. These exponents can be 
calculated analytically, by a suitable generalization of a technique devised 
by Ernst and Brito [29] . The power law decay of cooling mixtures has been 
shown to vary between 2 and oo. 

4 The one-dimensional gas 

The solution of the IMM given in (11) has nothing to do with the real dynam- 
ics of a cooling granular gas in the homogeneous regime, where no algebraic 
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Fig. 3. Asymptotic velocity distributions P(|«|,r) versus \v\/vo(t) for r = from 
the simulation of the Inelastic Maxwell Model in 2D (left) and 3D (right). In 2D 
the asymptotic distribution tail converges to the exactly computed exponent (model 
IMM-A). In 3D the convergence is slower, but the difference of the predicted ex- 
ponent for model IMM-A and IMM-B can be appreciated. In the bottom panels 
the convergences of the local logarithmic slope is put in evidence. (For the precise 
definition of the two models see [29].) Data refers to more than N = 10 6 particles. 



tails appear. Instead, a good agreement (in the HCS) is found between the so- 
lution of Boltzmann equation [36] and Molecular Dynamics [37] . However we 
will show in the following that the Inelastic Maxwell Model with the addiction 
of a special topology (the Id lattice) is able to reproduce the homogeneous 
regime and some essential features of the first instability (enhancement of 
velocity fluctuations). 



4.1 Molecular Dynamics 

An MD simulation of N hard rods on a ring of length L shows that, besides 
hydrodynamic instabilities discussed above, also a kinetic singularity appears 
that prevents long simulation runs, the so-called "inelastic collapse" [4,38]. 
The inelastic collapse is the divergence of the collision rate of a small set of 
grains, analogous to the divergence of the collision rate of a ball left bounc- 
ing on the floor. This means that a group of particles will experience an 
infinite number of collisions in a finite time. To avoid such a problem, many 
authors [9,39] have proposed the introduction of a velocity cut-off 5 such 
that, whenever two particles collide with a relative velocity smaller, in ab- 
solute value, than 6, the collision is elastic. Taking S 2 <C T min , being T min 
the minimum granular temperature that one expects to observe, one can be 
sure that the choice of 8 will not influence the main results of the simulation. 
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Fig. 4. Time behavior of kinetic energy for the Inelastic Hard Rods (left) and 
for the Inelastic Lattice Maxwell Model (right). The homogeneous Haff stage is 
evident only for quasi-elastic system, whereas the more inelastic ones enter almost 
immediately in the correlated regime. Note the different time units used (t is the 
physical time, defined only in MD simulations, while r is the cumulated number 
of collisions per particle). The Haff law E ~ t~ 2 ~ exp(— 7ot) is verified for both 
systems. The correlated regime presents a behavior t~ 2 ^ 3 independent of r for the 
Hard Rods, while appears diffusive (in collision units) r -1 ^ 2 and r-dependent for 
the Lattice Model. Data refer to N = 10 6 particles (both models). 

Performing an efficient Event Driven simulation scheme (based on a tree or- 
dering of the event times) it is possible to study easily very large systems, 
with N ~ 10 6 for thousands of collisions per particle. 

In the left panel of figure 4 the decay of kinetic energy E{t) = T{t) = {v 2 ) 
is shown. It displays, for different values of r, the same behavior with two 
regimes: in the HCS E(t) <~ t~ 2 (Haff's law), then it decays differently, i.e. 
E(t) <~ t~ 3 / 2 . Noticeably all the curves seem to collapse in the second regime, 
i.e. the asymptotic regime is equivalent (not only in the exponent, but also 
in the coefficients) for every value of r. 

The velocity pdf 's are shown in figure 5: in the HCS it appears the peculiar 
"two-peaks" form, which is also predicted by the study of the Boltzmann 
equation in the quasi-elastic limit [36] . In the second regime (which cannot be 
described by the Boltzmann equation, because of the presence of correlations) 
the pdf becomes a Gaussian. 

What can be said about correlations in the second regime? An inspec- 
tion of the velocity profile Vi vs. Xi (see top panel of figure 6) suggests the 
presence of shocks, i.e. asymmetric "discontinuities" of the velocity profile in 
correspondence with peaks of density (clusters) . A more quantitative informa- 
tion about the velocity field is given by the study of the structure function 
S(k,t) — v(k,t)v(—k,t) where v(k,t) is the Fourier transform of the field 
v(i,t), presented in the left frame of figure 7. A good collapse is obtained if 
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Homogeneous regime Asymptotic regime 




Fig. 5. Rescaled velocity distributions for the ID MD and in the ID ILMM, during 
the homogeneous (left) and the inhomogeneous phase (right). In the left frame is 
also shown the initial distribution (both models). The distributions refer to systems 
having the same energy. Data refer to N = 10 6 particles (both models) with r — 0.99 
and r = 0.5 (for the lattice model in the inhomogeneous regime). 
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Fig. 6. Portions of the instantaneous velocity profiles for the ID MD (top v(x,t), 
middle v(i,i)) and for the ID lattice gas model (bottom, v(i,r)). In the middle 
frame we display the MD profile against the particle label in order to compare the 
shocks and preshocks structures with the lattice gas model (the dotted lines show 
how shocks and preshocks transform in the two representations for the MD). Data 
refers to N — 2 ■ 10 4 particles, r = 0.99 (both models). 
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Fig. 7. Structure factors S(k, t) against fct 2//3 for the ID MD and against fcr 1 ^ 2 for 
the ID lattice gas model, in the inhomogeneous phase. Times are chosen so that the 
two systems have the same energies. Data refers to system with more than N = 10 
particles, r = 0.5 (both models) 



S(k,t) is plotted at different times t versus fci 2 / 3 - k/E{t). The presence of 
short scale defects (shocks) is put in evidence by the form of structure factor 
at high values of k, where S(k,t) ~ fc~ 2 : this tail of the structure function 
is expected, in d = 1, in the coarsening of a two-phases system, due to the 
non-analyticity (at short distances) of the spatial correlation function. This 
is also known as Porod law in the theory of phase ordering kinetics [40-42] . 

4.2 The Inelastic Lattice Maxwell Model 

The Inelastic Lattice Maxwell Model in d — 1 is defined as follows. A set of 
N scalar velocities are placed on N sites on a ring. After a time At = 2/N a 
couple of neighbors is chosen and is updated with the usual inelastic collision 
rule. The only constraint in the choice of the couple is that Uj > v i+ i\ this is 
called the kinematic constraint and represents the physical condition neces- 
sary for a collision. It must be noted that the choice of colliding particles does 
not depend on the modulus of the relative velocity: this is why we consider 
this a "Inelastic Maxwell Model" embedded in a Id lattice. In this model 
the only measure of time is r, i.e. the cumulated number of collisions per 
particle. This makes many measures difficult to be compared with analogous 
MD measurements: in fact in MD simulations r(t) depends upon the choice 
of 6 (the elastic cut-off), while all observables depend on t independently of 
5. 

The analysis of the ILMM follows closely that of MD simulations and the 
results [13,14] are presented, in parallel, in the same figures 4, 5, 6, 7. 
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The study of the decay of energy (figure 4) indicates that the model well 
reproduces the Haff regime (E(t) ~ cxp(— 7ot)) and displays a decay in the 
second regime of the form ~ t~ 1//2 . 

The velocity pdf 's reproduce exactly the ones measured in MD simulations 
(see fig. 5), showing the "two-peaks" form in the Haff stage and the Gaussian 
in the late stage. It must be stressed again that these features are obtained 
simply applying to the Inelastic Maxwell Model (which displays tails P(v) ~ 
v~ 4 ) a two-neighbors (instead of A^-neighbors) topology. 

The inspection of velocity profiles must be carefully carried: in this model 
the particles cannot move and therefore a comparison with MD should be 
made considering Vi vs. i instead of X{. The good agreement between MD 
and ILMM can be appreciated in the two bottom panels of figure 6: with this 
choice of abscissa, the shocks observed in MD simulations appear reversed 
and smoothed, very similar to the ones observed in the ILMM. Analogously 
the structure factor presents similar features with MD, mainly a k~ 2 behavior 
at high values of k (short scales) which is a signature of topological defects. 
In the lattice case the collapse can be obtained plotting S(k, r) against fcr 1 / 2 , 
i.e. ~ k/E(r) as in the MD case. We have also verified that, removing the 
kinematic constraint, this structure disappears and the model behaves as 
governed by a simple diffusion equation. 

5 The two-dimensional gas 

The number of interesting and realistic results obtained with the one di- 
mensional Inelastic Lattice Maxwell Model, in spite of its extreme simplicity, 
convinced us to extend our study to two dimensions. 

However, a smaller number of results is known from Molecular Dynamics 
in d = 2, especially in the correlated regime, which requires very long sim- 
ulations. Therefore the comparison between MD and the ILMM to the case 
d = 2 is more difficult. Only recently MD simulations start to give reliable 
measurements even for long times, and, where possible, the comparison with 
our ILMM is promising. 

More speculatively, we tried a series of less standard measurements with 
two different aims: the first is to try to understand what kind of effective 
equation could describe the velocity field of the ILMM; the second is to 
discuss some basic hypothesis necessary to a hydrodynamic description of 
this granular fluid. After a brief review of the known results from MD, we 
present the results obtained with our 2D ILMM. First the same analysis of 
energy decay, velocity pdf and structure factors is carried on, as in the Id case. 
Moreover, to better capture the nature of the topological defects generated 
by the dynamics, we consider the pdf of the velocity increments. Then we 
present an analysis of temporal correlations, mainly aimed to identify an 
effective equation for the velocity field. Finally we suggest a measurement 
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of a hydrodynamic temperature and we sketch a discussion about the very 
presence of scale separation. 

5.1 Known results from Molecular Dynamics 

To our knowledge few MD simulations have inspected the behavior of in- 
elastic hard disks in the late stage, i.e. after the end of Haff regime and the 
growth of instabilities. The departure from the homogeneity has been studied 
by Ernst and co-workers [6-8] and their MD simulations have shown results 
compatible with their predicted t 1 / 2 scaling law for the linear size of struc- 
tures (vortices and clusters). A figure in a work of Huthmann et al. [43] was 
the first suggestion that the global velocity pdf asymptotically returns to a 
Gaussian (while it has fatter tails in the HCS) and recently different works 
on large MD simulations have confirmed this [44,10]. Recently this "return- 
to-the-Gaussian" scenario was put in strict relation with the presence of a 
Burgers-like dynamics [10]. The debate on the asymptotic energy seems to 
be still open, even if there are strong indications of a universal (independent 
of r) decay of the form E(t) <~ i -1 [10]. A recent work of Isobe [20] demon- 
strates that many features of the "last" stage resemble the dynamics of 2D 
turbulence. 

5.2 The Inelastic Maxwell Lattice Model 

The ILMM in d — 2 is identical to that in d = 1 except for the fact that the 
particles have a 2-components vector velocity Vj and that they are placed on a 
triangular 2D lattice, i.e. every particle has 6 neighbors. The dynamics is the 
same as in d — 1; the kinematic constraint here reads (vj — Vj) • (r< — r,-) < 0. 

The decay of energy E{t) is shown in figure 8 for two different values of 
r. The Haff regime E(t) <~ exp(— 7 t) is well reproduced (when r = 0.2 it 
terminates too early to be appreciated). The asymptotic decay reads E{t) ~ 
r _1 . Both this measurement and the decay observed in the d = 1 case, are 
compatible with a universal law E(t) ~ r _d / 2 as for a purely diffusive field. 

The rescaled velocity pdf's (see figure 9), which are initially Gaussian, 
become non-Gaussian in the Haff regime, with larger tails, but return Gaus- 
sian in the late stage. The analysis that follows shows that, in the late stage, 
the velocity field is far from being homogeneous: the global velocity pdf is 
strongly influenced by this homogeneity and its Gaussian form could be sim- 
ply the signature of the presence of many independent domains (i.e. it is just 
the distribution of the average velocities of these domains). 

In the figures 10 and 11 we show the presence of vortices in the velocity 
fields. The second figure represents the same system at a later time, demon- 
strating the coarsening of vortices (their number reduces and their size in- 
creases) . 

In d = 2 the structure function is defined as 



16 Andrea Baldassarri et al. 




Fig. 8. Energy decay for r = 0.9 and r = 0.2 (1024 2 sites). The bold dashed line 
~ 1/t is a guide to the eye for the asymptotic energy decay, while the light dashed 
line is the exponential fit corresponding to the Haff law exp(— 2^/or). The Haff 
regime is too short to be observed in the system with r = 0.2. The two indicated 
times r = 24 and r = 860 correspond to the plots of fig. 16. 
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Fig. 9. Distributions of horizontal velocity at different times for the same system 
with r = 0.2 and N = 512 x 512. The distributions are rescaled in order to have unit 
variance. The initial distribution is a Gaussian. The distribution becomes broader 
in the uncorrelated phase (first regime), and then turns back toward a Gaussian. 
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Fig. 10. A (zoomed) snapshot of the velocity field at time r = 52 for the Inelastic 
Lattice Gas, d — 2, with r — 0.7 and size N = 512 x 512. The time has been chosen 
at the beginning of the correlated regime. It is evident the presence of vortices. All 
the velocities have been rescaled to arbitrary units, in order to be visible. 




Fig. 11. Another snapshot of the same system of fig. 10 but at a later time r = 535. 
The diameter of the vortices has grown up. All the velocities have been rescaled to 
arbitrary units, in order to be visible. 



^(fc,t)-^v t ' 1 (k,t)v^(-k,i) (12) 

k 

where the superscripts t, I indicate the transverse and longitudinal compo- 
nents of the field with respect to the wave vector k and the sum is over 
a circular shell of radius k. The analysis of the structure factor, shown in 
figure 12, is performed for two systems at different inelasticities. The com- 
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Fig. 12. Data collapse of the transverse (£') and longitudinal (S l ) structure func- 
tions for r = 0. and r = 0.9 (system size 1024 2 sites, times ranging from r = 500 to 
r = 10 4 ). The wave-number k has been multiplied by i/t. Notice the presence of 
the plateaux for the more elastic system. For comparison we have drawn the laws 
x~ 4 and exp(— x 2 ). 



mon feature is a good collapse of different time curves if plotted against 
fcr 1 / 2 which is the signature of a domain growth (diameter of the vortices) 
L(t) ~ t 1 / 2 . The elasticity however changes dramatically the other features: 
the more elastic system has a Gaussian structure factor at large scale (low 
values of k) and a plateau at small scales, i.e. quasi-elastic collisions random- 
ize neighbour velocities while the effect of inelasticity induces a structure 
at large scale. The less elastic system does not present such a plateau, but 
a decay fc -4 at small scales which is analogous to the Porod law [40,41] 
expected from the presence of defects in a phase ordering process. The differ- 
ence between transversal and longitudinal component (more appreciable in 
the quasi-elastic system) is coherent with the theoretical analysis discussed 
in the introduction, which predicts an earlier and stronger correlation for the 
transversal compoenent. 

In the d — 2 ILMM we have also observed shock-like phenomena. The 
inspection of the distribution of longitudinal and transversal velocity differ- 
ences, defined as 

A(R,i) = (v i+il -Vi)-^ (13a) 
4 t (R,i) = (v <+Ji -v<)x~ (13b) 
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Fig. 13. Probability densities of the longitudinal and transverse velocity incre- 
ments. The main figure shows the p.d.f. of the velocity differences for R = 1. The 
inset shows the Gaussian shape measured for R — 40 (larger than L(t) for this 
simulation: r = 0.2, t — 620, system size 2048 2 ). 



is contained in figure 13. The main graph shows the case R = 1, i.e. the distri- 
butions of gradients. Both the longitudinal and transverval components are 
distributed with non-Gaussian tails, but the longitudinal component presents 
also a strong asymmetry: this means that preferentially there are strong ve- 
locity differences in the direction of growing longitudinal velocity. This is 
compatible with the analysis of shocks in the bottom panels of figure 6 in the 
d = 1 case. When R^> 1 the distribution of velocity differences is a Gaussian. 

To better study the dynamics of the system we have calculated the two- 
times self-correlation of the velocity components 

C(T U T 2 ) = ^pM (14 ) 

During a short time transient, the self-correlation function of our model de- 
pends on n — T2, i.e. it is time translational invariant (TTI). Later, C(t\, t%) 
reaches an "aging" regime and depends only on the ratio x — t\/t%. This 
TTI transient regime is similar to what occurs during the coarsening process 
of a quenched magnetic system: the self-correlation of the local magnetiza- 
tion a{r w , t w + t) for r -C t w shows a TTI decay toward a constant value 
m 1.q{T quench) that is the square of the equilibrium magnetization, indicating 
that the local magnetization is evolving in an ergodic-like fashion. Thereafter 
the self-correlation decays with the aging scaling law indicated above. Obvi- 
ously, when T quenc h — > 0, the TTI transient regime disappears. In our model 
the behavior of the self-correlation is even more subtle, as the cooling pro- 
cess imposes a (slowly) decreasing "equilibrium" temperature T quenc h — > 0: 
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Fig. 14. Angular auto-correlation function A(t, t w ) for different values of the wait- 
ing time t w and r — 0.9 (f024 2 sites). The graph on the left shows the convergence 
to the t/t w diffusive scaling regime, for large t w . For small t w , a local minimum is 
visible (for such a quasi elastic dynamics). In the graph on the right the same data 
are plotted vst-t w : note that the small t w curves tend to collapse. For higher t w 
the position of the local minimum does not move sensibly, but its value grows and 
goes to 1 for large t w 

this progressively erodes the TTI regime and better resembles a finite rate 
quench. The same dependence on the TTI manifests itself in the angular 
auto-correlation, shown in figure 14: 



The non-monotonic behavior of A(t, t w ) suggests that the initial direction 
of the velocity induces a change in the velocities of the surrounding particles, 
which in turn generates, through a sequence of correlated collisions, a kind of 
retarded field oriented as the initial velocity. As t w increases the maximum 
is less and less pronounced. 

The behavior of the energy decay (E(r) ~ T~ d / 2 ), the law of growth of do- 
mains (L(t) ~ r 1 / 2 ) and the "aging" form of the two-times self-correlations, 
as well as the Gaussian velocity pdf 's and the Gaussian form at large scale 
in the structure factors, are strong evidence that the model behaves similarly 
to a system governed by the diffusion equation: 



A{t,t w ) = —y^cos(6j(T w +r) - 8i(T w )). 



(15) 




(16) 
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Fig. 15. Measure of the persistence of the Inelastic Lattice Gas in d = 2, for two 
different values of the restitution coefficient. The power-law regime corresponds to 
the correlated regime. The measured exponent 9 = 0.18 is, with this precision, equal 
to the persistence exponent numerically measured in the diffusive dynamics. 

where </> is a suitable field depending on the velocity field (it could be a scalar, 
for example the z-component of the vorticity field V x v) . 

To perform a further check of this conjecture, we have measured the per- 
sistence exponent [45] of the model, which is very sensitive to the class of 
universality. We have counted the number N s (r) of sites where the x velocity 
component never changed sign from the starting time of the dynamics up 
to time r. In the correlated regime we observe N s (t) <~ t~ with 9 = 0.18 
which is equal, up to this precision, to the persistence exponent numerically 
measured for the diffusive dynamics. 

In the spirit of verifying the possibility of a hydrodynamic (mesoscopic) 
description of this model, we have also analysed the average local granular 
temperature T ai defined as: 

T ff = (|v-(v) ff | 2 ) ff (17) 

where (■■■)„ means an average on a region of linear size a. 

If we call L{t) a characteristic correlation length of the system, since 
when a 3> L(t) the local average tends to the global (zero) momentum, then 
linv-.oc T a = E. For a < L{t), instead, T a < E. The behavior of T a in 
the uncorrelated (Haff) regime and in the correlated (asymptotic) regime for 
two different values of r is presented in figure 16. For quasi elastic systems 
T a exhibits a plateau for 1 <C o <C L{t) that identifies the strength of the 
internal noise (see also the plateau in the structure factor, figure 12) and 
indicates the mesoscopic scale necessary for a hydrodynamics description. 
The local temperature ceases to be well defined for smaller r, revealing an 
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Fig. 16. The scale dependent temperature, T a , defined as function of the coarse 
graining size a for r = 24 (incoherent regime) and r = 860 (correlated regime for 
both choices of r): the system is the same of figure 8. In the correlated regime 
the more elastic case presents a plateau at intermediate wavelengths, indicating a 
well defined mesoscopic temperature, and therefore a clear separation between the 
microscopic and the macroscopic scales. 

absence of scale separation between microscopic and macroscopic fluctuations 
in the strongly inelastic regime [46]. 

6 Conclusions 

We have analysed essentially three different problems with a strong intercon- 
nection, in the framework of the physics of cooling granular gases. Our first 
step (IMM) was to demonstrate that a simplified ID version of the Boltzmann 
equation for inelastic gases can be exactly solved by a P(v) with power law 
tails. However, realistic granular systems never show such large tails: the key 
feature that must be added to this oversimplified model is a realistic topol- 
ogy. Therefore our second step (ILMM, d = 1) was to embed the IMM on 
a ID lattice. The P(v) changes dramatically and becomes identical to those 
measured in MD simulations, also in the velocity-correlated stage. Moreover, 
the ILMM allows for a study of spatial correlations, resulting in shock-like 
structures very similar to the one observed in MD and strong analogies with 
the MD structure functions. Our last step (ILMM, d — 2) has been the 
passage from d = 1 to d = 2, which mainly results in the appearance of vor- 
tices. Vortices are ubiquitous in granular gases experiments and simulations. 
The ILMM simply produces vortices as a result of the competion between 
the parallelizing effect of inelastic collisions and the constraint given by the 
conservation of total momentum. Moreover, a large set of evidences (energy 
decay, structure factors, Gaussian velocity pdf's, aging self-correlations, per- 
sistence exponent) indicate that the dynamical behavior of the system is 
compatible with a diffusion equation, even if short scale defects (appearing 
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as shocks, internal noise and tails a la Porod in the structure functions) make 
this model richer. 
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